library(stats)
suppressPackageStartupMessages(
  library(fpp2)
)
library(forecast)
suppressPackageStartupMessages(
  library(dynlm)
)
library(magrittr)
suppressPackageStartupMessages(
  library(tseries)
)
L=read.csv('LaborForce.csv',header=T)

Labor = ts(L$LNU01327662,start=1992,frequency=12)

autoplot(Labor,main="US Labor Force Participation Rate",ylab="Percent",xlab="Year")

acf(Labor)

pacf(Labor)

t = time(Labor)
model1=dynlm(Labor~t)
plot1=ts(predict(model1),start=1992,frequency=12)
autoplot(plot1)+autolayer(Labor)

model2=dynlm(Labor~t+sin(2*pi*t/13)+cos(2*pi*t/13))
plot2=ts(predict(model2),start=1992,frequency=12)
autoplot(plot2) + autolayer(Labor)

model3=dynlm(Labor~t+I(t*cos(2*pi*t/15))+I(t*sin(2*pi*t/15)))
plot3=ts(predict(model3),start=1992,frequency=12)
autoplot(plot3) + autolayer(Labor)

model4=dynlm(Labor~t+sin(2*pi*t/13)+t*cos(2*pi*t/13))
plot4=ts(predict(model4),start=1992,frequency=12)
autoplot(plot4) + autolayer(Labor)

model5=dynlm(Labor~t+sin(2*pi*t/14)+t*cos(2*pi*t/14))
plot5=ts(predict(model5),start=1992,frequency=12)
autoplot(plot5) + autolayer(Labor)

model6=dynlm(Labor~t+cos(2*pi*t/12))
plot6=ts(predict(model6),start=1992,frequency=12)
autoplot(plot6) + autolayer(Labor)

model7=dynlm(Labor~t+t*cos(2*pi*t/15)+t*sin(2*pi*t/15))
plot7=ts(predict(model7),start=1992,frequency=12)
autoplot(plot7) + autolayer(Labor)

model8=dynlm(Labor~I(t*t)+t*cos(2*pi*t/13)+sin(2*pi*t/13))
plot8=ts(predict(model8),start=1992,frequency=12)
autoplot(plot8) + autolayer(Labor)

t_adj=t-1992
model28=nls(Labor~a+b*t_adj+
             c*cos(2*pi*t_adj/p)+d*sin(2*pi*t_adj/p)+
             e*t_adj*cos(2*pi*t_adj/p)+f*t_adj*sin(2*pi*t_adj/p)+
             g*cos(4*pi*t_adj/p)+h*sin(4*pi*t_adj/p)+
             i*t_adj*cos(4*pi*t_adj/p)+j*t_adj*sin(4*pi*t_adj/p),
           start=data.frame(a=82,b=-0.29,p=14,
                            c=-0.0934,d=0.28,
                            e=0.031317,f=-0.035067,
                            g=0.01,h=0.01,
                            i=0.01,j=0.01))
plot28=ts(predict(model28),start=1992,frequency=12)
autoplot(plot28) + autolayer(Labor)

t_adj=t-1992
model29=nls(Labor~a+b*t_adj+
             c*cos(2*pi*t_adj/p)+d*sin(2*pi*t_adj/p)+
             e*t_adj*cos(2*pi*t_adj/p)+f*t_adj*sin(2*pi*t_adj/p)+
             g*cos(4*pi*t_adj/p)+h*sin(4*pi*t_adj/p),
           start=data.frame(a=82,b=-0.29,p=14,
                            c=-0.0934,d=0.28,
                            e=0.031317,f=-0.035067,
                            g=0.01,h=0.01))
plot29=ts(predict(model29),start=1992,frequency=12)
autoplot(plot29) + autolayer(Labor)

t_adj=t-1992
model30=nls(Labor~a+b*t_adj+
             c*cos(2*pi*t_adj/p)+
             e*t_adj*cos(2*pi*t_adj/p)+
             g*cos(4*pi*t_adj/p),
           start=data.frame(a=82,b=-0.29,p=14,
                            c=-0.0934,
                            e=0.031317,
                            g=0.01))
plot30=ts(predict(model30),start=1992,frequency=12)
autoplot(plot30) + autolayer(Labor)

model31=dynlm(Labor~t+
                t*cos(2*pi*t/15.3)+t*sin(2*pi*t/15.3)+
                cos(4*pi*t/15.3)+sin(4*pi*t/15.3))
plot31=ts(predict(model31),start=1992,frequency=12)
autoplot(plot31) + autolayer(Labor)

summary(model1)
## 
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
## 
## Call:
## dynlm(formula = Labor ~ t)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.38458 -0.35754  0.00369  0.36341  1.41844 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 680.603408   7.270444   93.61   <2e-16 ***
## t            -0.300534   0.003624  -82.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.537 on 334 degrees of freedom
## Multiple R-squared:  0.9537, Adjusted R-squared:  0.9535 
## F-statistic:  6876 on 1 and 334 DF,  p-value: < 2.2e-16
summary(model2)
## 
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
## 
## Call:
## dynlm(formula = Labor ~ t + sin(2 * pi * t/13) + cos(2 * pi * 
##     t/13))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.03689 -0.29964  0.00565  0.33782  1.10905 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        657.877547   6.356894 103.490  < 2e-16 ***
## t                   -0.289201   0.003169 -91.260  < 2e-16 ***
## sin(2 * pi * t/13)  -0.258539   0.034085  -7.585 3.37e-13 ***
## cos(2 * pi * t/13)  -0.372962   0.036977 -10.086  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4471 on 332 degrees of freedom
## Multiple R-squared:  0.9681, Adjusted R-squared:  0.9678 
## F-statistic:  3356 on 3 and 332 DF,  p-value: < 2.2e-16
summary(model3)
## 
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
## 
## Call:
## dynlm(formula = Labor ~ t + I(t * cos(2 * pi * t/15)) + I(t * 
##     sin(2 * pi * t/15)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1393 -0.3780 -0.0018  0.3711  1.1479 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                6.553e+02  7.144e+00  91.725   <2e-16 ***
## t                         -2.879e-01  3.562e-03 -80.843   <2e-16 ***
## I(t * cos(2 * pi * t/15))  1.738e-04  1.967e-05   8.835   <2e-16 ***
## I(t * sin(2 * pi * t/15)) -2.015e-05  1.934e-05  -1.042    0.298    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4841 on 332 degrees of freedom
## Multiple R-squared:  0.9626, Adjusted R-squared:  0.9622 
## F-statistic:  2847 on 3 and 332 DF,  p-value: < 2.2e-16
summary(model4)
## 
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
## 
## Call:
## dynlm(formula = Labor ~ t + sin(2 * pi * t/13) + t * cos(2 * 
##     pi * t/13))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.08380 -0.29193  0.02731  0.29376  0.97821 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          664.145333   6.140604 108.156  < 2e-16 ***
## t                     -0.292303   0.003061 -95.507  < 2e-16 ***
## sin(2 * pi * t/13)    -0.249624   0.032473  -7.687 1.73e-13 ***
## cos(2 * pi * t/13)    55.931373   9.444133   5.922 7.95e-09 ***
## t:cos(2 * pi * t/13)  -0.028067   0.004708  -5.962 6.39e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4255 on 331 degrees of freedom
## Multiple R-squared:  0.9712, Adjusted R-squared:  0.9708 
## F-statistic:  2787 on 4 and 331 DF,  p-value: < 2.2e-16
summary(model5)
## 
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
## 
## Call:
## dynlm(formula = Labor ~ t + sin(2 * pi * t/14) + t * cos(2 * 
##     pi * t/14))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.14586 -0.30482  0.02367  0.29644  0.90439 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          654.846505   6.154874 106.395  < 2e-16 ***
## t                     -0.287650   0.003068 -93.756  < 2e-16 ***
## sin(2 * pi * t/14)    -0.107086   0.032726  -3.272  0.00118 ** 
## cos(2 * pi * t/14)    79.061920   8.956710   8.827  < 2e-16 ***
## t:cos(2 * pi * t/14)  -0.039613   0.004466  -8.870  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4183 on 331 degrees of freedom
## Multiple R-squared:  0.9721, Adjusted R-squared:  0.9718 
## F-statistic:  2888 on 4 and 331 DF,  p-value: < 2.2e-16
summary(model6)
## 
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
## 
## Call:
## dynlm(formula = Labor ~ t + cos(2 * pi * t/12))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.15791 -0.32242  0.02468  0.29195  1.14316 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        672.40224    6.09863  110.25   <2e-16 ***
## t                   -0.29643    0.00304  -97.50   <2e-16 ***
## cos(2 * pi * t/12)  -0.42973    0.03537  -12.15   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4477 on 333 degrees of freedom
## Multiple R-squared:  0.9679, Adjusted R-squared:  0.9677 
## F-statistic:  5020 on 2 and 333 DF,  p-value: < 2.2e-16
summary(model7)
## 
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
## 
## Call:
## dynlm(formula = Labor ~ t + t * cos(2 * pi * t/15) + t * sin(2 * 
##     pi * t/15))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.14894 -0.26691  0.05073  0.30602  0.90700 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          650.602081   6.119585 106.315  < 2e-16 ***
## t                     -0.285535   0.003051 -93.597  < 2e-16 ***
## cos(2 * pi * t/15)   -62.443236   8.039196  -7.767 1.02e-13 ***
## sin(2 * pi * t/15)    70.253341   8.698830   8.076 1.27e-14 ***
## t:cos(2 * pi * t/15)   0.031317   0.004007   7.816 7.39e-14 ***
## t:sin(2 * pi * t/15)  -0.035067   0.004337  -8.085 1.19e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.412 on 330 degrees of freedom
## Multiple R-squared:  0.9731, Adjusted R-squared:  0.9726 
## F-statistic:  2383 on 5 and 330 DF,  p-value: < 2.2e-16
summary(model8)
## 
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
## 
## Call:
## dynlm(formula = Labor ~ I(t * t) + t * cos(2 * pi * t/13) + sin(2 * 
##     pi * t/13))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.13848 -0.27107  0.05675  0.29772  0.91265 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -8.797e+03  1.826e+03  -4.816 2.23e-06 ***
## I(t * t)             -2.351e-03  4.538e-04  -5.180 3.87e-07 ***
## t                     9.139e+00  1.821e+00   5.019 8.49e-07 ***
## cos(2 * pi * t/13)    3.448e+01  9.994e+00   3.451 0.000632 ***
## sin(2 * pi * t/13)   -1.858e-01  3.362e-02  -5.526 6.67e-08 ***
## t:cos(2 * pi * t/13) -1.739e-02  4.981e-03  -3.491 0.000547 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4098 on 330 degrees of freedom
## Multiple R-squared:  0.9733, Adjusted R-squared:  0.9729 
## F-statistic:  2409 on 5 and 330 DF,  p-value: < 2.2e-16
summary(model28)
## 
## Formula: Labor ~ a + b * t_adj + c * cos(2 * pi * t_adj/p) + d * sin(2 * 
##     pi * t_adj/p) + e * t_adj * cos(2 * pi * t_adj/p) + f * t_adj * 
##     sin(2 * pi * t_adj/p) + g * cos(4 * pi * t_adj/p) + h * sin(4 * 
##     pi * t_adj/p) + i * t_adj * cos(4 * pi * t_adj/p) + j * t_adj * 
##     sin(4 * pi * t_adj/p)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## a 81.721164   0.284669 287.074  < 2e-16 ***
## b -0.278969   0.018513 -15.069  < 2e-16 ***
## p 15.624579   2.453688   6.368 6.54e-10 ***
## c -0.509344   0.122191  -4.168 3.94e-05 ***
## d  0.297235   0.730436   0.407   0.6843    
## e  0.057260   0.022125   2.588   0.0101 *  
## f -0.003323   0.081301  -0.041   0.9674    
## g  0.154910   0.241634   0.641   0.5219    
## h  0.054556   0.084684   0.644   0.5199    
## i -0.001288   0.014332  -0.090   0.9285    
## j -0.002584   0.013551  -0.191   0.8489    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.405 on 325 degrees of freedom
## 
## Number of iterations to convergence: 17 
## Achieved convergence tolerance: 6.825e-06
summary(model29)
## 
## Formula: Labor ~ a + b * t_adj + c * cos(2 * pi * t_adj/p) + d * sin(2 * 
##     pi * t_adj/p) + e * t_adj * cos(2 * pi * t_adj/p) + f * t_adj * 
##     sin(2 * pi * t_adj/p) + g * cos(4 * pi * t_adj/p) + h * sin(4 * 
##     pi * t_adj/p)
## 
## Parameters:
##    Estimate Std. Error  t value Pr(>|t|)    
## a 81.756853   0.063446 1288.595  < 2e-16 ***
## b -0.281288   0.004088  -68.812  < 2e-16 ***
## p 15.305461   0.525971   29.099  < 2e-16 ***
## c -0.488013   0.069050   -7.067 9.58e-12 ***
## d  0.205529   0.149410    1.376 0.169886    
## e  0.053475   0.007593    7.042 1.12e-11 ***
## f  0.006954   0.016004    0.435 0.664199    
## g  0.127955   0.038381    3.334 0.000955 ***
## h  0.045083   0.056791    0.794 0.427860    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4039 on 327 degrees of freedom
## 
## Number of iterations to convergence: 9 
## Achieved convergence tolerance: 3.75e-06
summary(model30)
## 
## Formula: Labor ~ a + b * t_adj + c * cos(2 * pi * t_adj/p) + e * t_adj * 
##     cos(2 * pi * t_adj/p) + g * cos(4 * pi * t_adj/p)
## 
## Parameters:
##    Estimate Std. Error  t value Pr(>|t|)    
## a 81.825556   0.047067 1738.479  < 2e-16 ***
## b -0.284602   0.003076  -92.517  < 2e-16 ***
## p 16.165036   0.132575  121.931  < 2e-16 ***
## c -0.402314   0.066260   -6.072 3.48e-09 ***
## e  0.051564   0.004396   11.730  < 2e-16 ***
## g  0.099924   0.033334    2.998  0.00293 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4236 on 330 degrees of freedom
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 3.962e-06
summary(model31)
## 
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
## 
## Call:
## dynlm(formula = Labor ~ t + t * cos(2 * pi * t/15.3) + t * sin(2 * 
##     pi * t/15.3) + cos(4 * pi * t/15.3) + sin(4 * pi * t/15.3))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.09435 -0.25024  0.04469  0.28124  0.95038 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             6.421e+02  6.201e+00 103.546  < 2e-16 ***
## t                      -2.813e-01  3.092e-03 -90.996  < 2e-16 ***
## cos(2 * pi * t/15.3)   -2.235e+01  8.887e+00  -2.515 0.012373 *  
## sin(2 * pi * t/15.3)   -1.054e+02  8.534e+00 -12.356  < 2e-16 ***
## cos(4 * pi * t/15.3)   -1.281e-01  3.309e-02  -3.871 0.000131 ***
## sin(4 * pi * t/15.3)    4.461e-02  3.320e-02   1.343 0.180046    
## t:cos(2 * pi * t/15.3)  1.104e-02  4.428e-03   2.494 0.013121 *  
## t:sin(2 * pi * t/15.3)  5.274e-02  4.257e-03  12.389  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4032 on 328 degrees of freedom
## Multiple R-squared:  0.9743, Adjusted R-squared:  0.9738 
## F-statistic:  1780 on 7 and 328 DF,  p-value: < 2.2e-16
AIC(model1,model2,model3,model4,model5,model6,model7,model8,
    model28,model29,model30,model31)
BIC(model1,model2,model3,model4,model5,model6,model7,model8,
    model28,model29,model30,model31)
checkresiduals(model1)

## 
##  Breusch-Godfrey test for serial correlation of order up to 24
## 
## data:  Residuals
## LM test = 266.34, df = 24, p-value < 2.2e-16
checkresiduals(model2)

## 
##  Breusch-Godfrey test for serial correlation of order up to 24
## 
## data:  Residuals
## LM test = 242.49, df = 24, p-value < 2.2e-16
checkresiduals(model3)

## 
##  Breusch-Godfrey test for serial correlation of order up to 24
## 
## data:  Residuals
## LM test = 251.85, df = 24, p-value < 2.2e-16
checkresiduals(model4)

## 
##  Breusch-Godfrey test for serial correlation of order up to 24
## 
## data:  Residuals
## LM test = 233.31, df = 24, p-value < 2.2e-16
checkresiduals(model5)

## 
##  Breusch-Godfrey test for serial correlation of order up to 24
## 
## data:  Residuals
## LM test = 230.37, df = 24, p-value < 2.2e-16
checkresiduals(model6)

## 
##  Breusch-Godfrey test for serial correlation of order up to 24
## 
## data:  Residuals
## LM test = 243.2, df = 24, p-value < 2.2e-16
checkresiduals(model7)

## 
##  Breusch-Godfrey test for serial correlation of order up to 24
## 
## data:  Residuals
## LM test = 229.23, df = 24, p-value < 2.2e-16
checkresiduals(model8)

## 
##  Breusch-Godfrey test for serial correlation of order up to 24
## 
## data:  Residuals
## LM test = 228.81, df = 24, p-value < 2.2e-16
checkresiduals(model28)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

checkresiduals(model29)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

checkresiduals(model30)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

checkresiduals(model31)

## 
##  Breusch-Godfrey test for serial correlation of order up to 24
## 
## data:  Residuals
## LM test = 232.38, df = 24, p-value < 2.2e-16